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ABSTRACT 


Despite si gnificant achievements in computational fluid dynamics, there still remain many fluid 
flow phenomena not well understood. For example, the prediction of temperature distributions is 
inamiratft when temperature gradients are high, particularly in shock wave turbulent boundary 
layer interactions close to the wall. Complexities of fluid flow phenomena include transition to 
turbulence, relaminarization, separated flows, transition between viscous and invisdd, 
incompressible and compressible flows, among others, in all speed regimes. The purpose of this 
paper is to introduce a new approach, called the Flowfield-Dependent Mixed Explicit-Implicit 
(FDMEI) method, in an attempt to resolve these difficult issues in CFD. In this process, a total of 
sue implicitness parameters characteristic of the current flowfield are introduced. They are 
from the current flowfield or changes of Mach numbers, Reynolds numbers, Peclet 
numb ers, and Damkohler numbers (if reacting) at each nodal point and time step. This implies that 
every nodal point or element is provided with different or unique numerical scheme according to 
their current flowfield situations, whether compressible, incompressible, viscous, inviscid, laminar, 
turbulent, reacting, or nonreacting. In this procedure, discontinuities or fluctuations of all 
variables between adjacent nodal points are determined accurately. If these implicitness 
parameters are fixed to certain numbers instead of being calculated from the flowfield information, 
fh#»n practically all currently available schemes of finite differences or finite elements arise as 
sp eci al cases. Some benchmark problems to be presented in this paper will show the validity, 
accuracy, and efficiency of the proposed methodology. 
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1. INTRODUCTION 

Nearly half a century has elapsed since the digital computer revolutionized computational 
technologies in engineering and mathematical physics. During this time finite difference methods 
(FDM) have dominated the field of computational fluid dynamics (CFD) [1-7], whereas the 
opposite is true for finite element methods (FEM) in solid mechanics. In recent years, however, 
the trend toward finite element methods in CFD appears to be increasingly favorable [8-14]. 

In general, the analyst preoccupied with the methods of his choice based on his 
educational background or research experience is seldom motivated to investigate other options. 
Thus, today the gap between these two disciplines is widely apart, despite the fact that the 
thorough understanding of the relations between FDM and FEM is beneficial. The purpose of this 
paper is an attempt to call for a new approach in which both FEM and FDM can be united toward 
the common goal of achieving the highest level of accuracy and efficiency in CFD. Similarities 
and dissi milarities must be identified in order to recognize merits and demerits of each method and 
to enable the analysts to choose the most desirable approach suitable for the particular task at «* 
hand. 


One of the most important questions in CFD is how to deal with large gradients of the 
variable (density, velocity, pressure, temperature, and source terms). Rapid changes of Mach 
numbers, Reynolds numbers, Peclet numbers, and Damkohler numbers (if reacting) between 
adjacent nodal points or elements can be a crucial factor in determining whether the chosen 
computational scheme will succeed or faiL Furthermore, proper treatments for incompressibility 
and compressibility, viscous and inviscid flows, subsonic and supersonic flows, laminar and 
turbulent flows, nonreacting and reacting flows are extremely important. The most general case 
of fluid dynamics where these various flow properties may be depicted in external and internal 
hypersonic flows is shown in Fig la,b. A typical reacting flow (hydrogen-air reaction) can also be 
seen in Fig. lc. 

Can a single formulation and computer code be made available to satisfy all the 
requirements mentioned above? Can a single mathematical formulation lead to most of the 
currently available computational schemes both in FDM and FEM as special cases? Most 
importantly, will such an approach guarantee accuracy and efficiency? In this paper, we respond 
to these questions positively, based on the results obtained through example problems. 

Toward this goal, our approach is based on the following procedure [15, 16], known as 
the Flowfield-Dependent Mixed Explicit-Implicit (FDMEI) scheme: 

(a) Write the Navier-Stokes system of equations in a conservation form. 

(b) Expand the conservation variable U" +1 in Taylor series up to and including the second- 
order time derivatives of the conservation variables. 

(c) Introduce in step (b) six different flowfield-dependent implicitness parameters which 
are calculated from the changes in Mach numbers, Reynolds numbers, Peclet numbers, 
and Damkohler numbers (if reacting) between nodal points or local elements. 
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(d) Substitute step (a) and (c) into step (b) to obtain the increments of the conservation 
variables ALT* 1 . As a result, the final form resembles the implicit factored scheme of 
Beam and Warming [1], but much more rigorous. 

(e) Step (d) may be used either in FDM or FEM. 

The computational procedure as described above is capable of resolving complex 
properties of fluid flows in general with shock waves, turbulence, and reacting flows in particular. 

(1) Incompressible flows are dependent on changes in Reynolds number between nodal 
points in FDM and within local elements in FEM. Incompressibility conditions are 
characterized by these changes in Reynolds number. 

(2) Compressible flows are dependent on changes in both Mach number and Reynolds 
n umb er between nodal points in FDM and within local elements in FEM. Dilatational 
dissipation is characterized by these changes in Mach number and Reynolds number. 

(3) Shock waves in compressible flows are dependent on changes in Mach number •* 
between nodal points in FDM and within local elements in FEM. Shock wave 
di sc ontinuities are characterized by these changes in Mach number. 

(4) High temperature gradient flows are dependent on changes in Peclet number between 
nodal points in FDM and within local elements in FEM. The convection vs diffusion 
in heat transfer is characterized by these changes in Peclet number. 

(5) Reacting flows are dependent on changes in Damkohler number between nodal points 
in FDM and within local elements in FEM. The mass source vs convective transfer, 
mass source vs diffusive transfer, heat source vs convective heat transfer, heat source 
vs conductive heat transfer, and heat source vs diffusive heat transfer are characterized 
by these changes in Damkohler number. 

(6) Direct numerical simulation (DNS1 for turbulent flows in which mesh refinements are 
carried out until turbulence length microscales are resolved without turbulence models 
can not be reliable particularly for high speed compressible turbulent flows unless the 
computational scheme is capable of treating high gradients of variables as described in 
(2) above. To improve turbulence calculations, Legendre polynomial spectral modes 
may be added as shown in [15]. Whether or not the spectral mode approach is 
advantageous for an overall computational efficiency remains to be seen. Due to the 
limita tion of computer time, the example problems in this paper are not intended for 
DNS microscale resolutions. 

Details of the mathematical formulations as described above are presented in Section 2, 
imp lemen tatio n and computational process in Section 3, some example problems in Section 4, 
and concluding remarks in Section 5. 
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2. MATHEMATICAL FORMULATIONS 


For the general purpose program considering the compressible viscous reacting flows, we 
write the conservation form of the Navier-Stokes system of equations as 


dt dx, dX; 


( 1 ) 


where U, F, G,, and B denote the conservation flow variables, convection flux variables, 
diffusion flux variables, and source terms, respectively. 
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where f- - YJt, is the body force, Y k is the chemical species, H\ is the zero-point enthalpy, 

w is the reaction rate, and is the binary diffiisivity. Additional equations for vibrational and 
electronic energies may be included in (1) for hypersonics. 

Expanding the conservation variables U in Taylor series including the first and second 
derivatives, we have 


U ~ = 0 . + ^ + ^^ +0( A,’) 

dt 2 dt 1 ' ' 

where s { and s 2 are the implicitness parameters defined such that 
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( 2 ) 

( 3 ) 

( 4 ) 


with Alf +1 = - U" . It is assumed that the convection flux F ( is a function of U and the 

diffusion fluxG ; is a function of both U and its gradient U ; . Thus, we have 


au _ 3F, 3G,- t ^ 

dt dx t dx, 


( 5 ) 

( 6 ) 


where the convection Jacobian a, , the diffusion Jacobian b ; , the diffusion gradient Jacobian c i; , 
and the source Jacobian d are defined as 
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a = 


aE 
au ’ 



c — £l , 
dU ; 



(7) 


Substituting (3) - (6) into (2) and assuming the product of the diffusion gradient Jacobian with 
third order spatial derivatives to be negligible, we obtain 
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( 8 ) 


In order to provide different implicitness (different numerical treatments or schemes) to 
different physical quantities, we reassign s 2 and s 2 associated with the diffusion and source terms, 
respectively. 


j,AG ; => r, AG ; 

, SjAB => s 5 AB 

(9a) 

r 2 AG, => r 4 AG ; 

, s 2 AB => s 6 AB 

(9b) 


with the various implicitness parameters defined as 

Sj = first order convection implicitness parameter 
s 2 = second order convection implicitness parameter 
s 3 = first order diffusion implicitness parameter 
s 4 = second order diffusion implicitness parameter 
s j = first order source term implicitness parameter 
s 6 = second order source term implicitness parameter 

The first order implicitness parameters s v s v and s 5 will be shown to be flowfield 
dependent with the solution accuracy assured by taking into account the flowfield gradients, 
whereas the second order implicitness parameters s 2 , s 4 , and s 6 , which are also flowfield 
dependent, mainly act as artificial viscosity, contributing to the solution stability. 

Substituting these implicitness parameters as defined in (9) into (8), we obtain 
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with 


9B 
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oU 


( 11 ) 


For simplicity we may rearrange (10) in a compact form. 
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Q ' = s: 
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We may allow the source terms in the LHS of (13) to lag from the time step n + 1 to n, so that 
(13b) can be written as . — 
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Note that the Beam- Warming Scheme [1] can be written in the form identical to (18) with 
the following definitions ofE ; , E~ and Q" 

E, = mAt (a ; + b ; ) , with m = 0/(1+^) (20) 

E^- = mAt c i; (21) 
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( 22 ) 


where the cross derivative terms appearing in Q" for the Beam- Warming scheme are included in 
the second derivative terms on the LHS. The Beam- Warming scheme is seen to be a special case 
of the FDMEI equations if we set s l = s 3 = m, s 2 = s 4 = s 5 = s 6 = 0, in (18), with adjustments of Q 
on the RHS. The stability analysis of the Beam-Warming scheme requires % > 0.385 and 0 = 2 + 
This will fix the implicitness parameter m to be 0.639 < m< 0.75. It can be shown that the 
FDMEI equations as derived in (9) are capable of producing practically all existing FDM and 
FEM schemes. Some examples are shown in Appendix A. 

Contrary to the Beam-Warming scheme, the FDMEI approach is to obtain the implicitness 
par am eters from the current flowfield variables at each and every nodal point rather than by fixing 
the implicitness parameters to certain predetermined numbers and using them for the entire flow 
domain irrespective of the local flowfield variation from one point to another. These implicitness 
parameters may be determined for spatial and temporal bases as depicted in Fig. 2. The final 
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values of implicitness parameters at any point and at any time can be obtained as the average of 
both spatial and temporal contributions: 

Convection Impli citness P arametersi 


s i 


rain(r,l) 
. 0 

1 


r > a 

r< a , *0 

=0 


with 


s 2 =s" , 0 < n < 1 (23) 


(24) 


where the maximum and minimum Mach numbers are calculated between adjacent nodal points in 
FDM or within a local finite element in FEM for spatial implicitness parameters (Fig. 2a) and 
between the time step at n and n - 1 for temporal implicitness parameters (Fig. 2b), and a is a 
us6r*sp£cificd sm all number (ot = 0.01). Here it is seen that Sj is directly related to the flowfield, 
whereas s 2 depends on s x such that s 2 = s " . The primary role of s x is to ensure the solution 
accuracy by properly accommodating the convection gradients, whereas that of s 2 is to act as 
artificial viscosity, for solution stability. 

Diffusion Implici tness Parameters: 


h 


rain(s,l) s>P 

0 s<p,Re mia ^0, orPe^O 

1 Re mia =0, orPe mi ,=0 


s 4 = s" , 0 < n < 1 


(25) 


with 

s = V Re^-Re^/Re^ or s = (26a,b) 


where the maximum and minimum Reynolds numbers or maximum and minimum Peclet numbers 
are calc ulated similarly as in s, for spatial and temporal implicitness parameters, and |3 is a user- 
specified small number (P = 0.01). If temperature gradients are large, it is possible that Peclet 
numbers instead of Reynolds numbers will dictate the diffusion implicitness parameters. The 
larger value of s 3 is to be chosen, as obtained either from (26a) or (26b). Note also that s 4 = s" 
with s } ensuring the solution accuracy by taking into account the diffusion gradients, and here 
again, s 4 plays the role of artificial viscosity, for solution stability. 

Source Term Implicit ness Parameters'. 

For the case of chemically reacting flows theDa (Damkohler number) must be used 
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frain(r,l) 
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s 6 =s; , 0< n < 1 (27) 


with 
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where the maximum and minimum Damkohler numbers are calculated similarly as in s l and s 3 for 
s patial and temporal implicitness parameters, and y is a user-specified small number (y = 0.01). 
The relationship between s 5 and s 6 is similar to those for convection and diffusion implicitness 
parameters such that s 6 = s” with s s and s 6 controlling the solution accuracy and solution stability, 
respectively. The average of both spatial and temporal implicitness parameters will be adopted for 
use in computations at any point (element) and time. 

Relationships between all physical phenomena and the corresponding numerical treatments 
are characterized by the balance between the first order implicitness parameters (s v s v s 5 ) and the 
second order implicitness parameters (s 2 , s 4 , s 6 ), ensuring the computational accuracy and 
computational stability, respectively. The idea is to provide adequate (no more and no less than 
required) amount of numerical viscosity in order to preserve the computational accuracy. Note 
that the defini tio ns for the second order implicitness parameters have been modified from those 
reported in [15, 16] in order to meet the above requirements (Fig. 3). Initially, it was thought that 
the second order implicitness parameters should be the direct opposite compliances of the first 
order implicitness parameters (s 2 = 1 - s v s 4 =1 - s v s 6 =1 - s 5 ) [16] such that the second order 
implicitness pa rame ters are the maximum and minimum, respectively, for the minimum and 
maximum values of the first order implicitness parameters. Unfortunately, such definition resulted 
in too little numerical viscosities for the high values of the first order implicitness parameters. 
Subsequently, the limiting values (0.5) of the second order implicitness parameters were provided 
such that s r = max(l - , 0.5), etc. as experimented in [15]. However, it was noted that both 

first and second order parameters should assume the same values at the both extremes at zero and 
unity with the second order implicitness parameters being reasonably large for all values of the 
first order implicitness parameters. Thus, the second order implicitness parameters given above 
are the nonlinear continuous functions of the first order implicitness parameters satisfying these 
requirements. The range of the constant n is 0 < n < 1, although rt = } has been found to be the 
optimum, exhibiting the best convergence rate for reasonably high CFL numbers in the example 
problems presented in Section 4. 

The flowfield dependent implicitness parameters as defined above are capable of allowing 
various numerical schemes to be automatically generated, as summarized below: 

(1) The first order implicitness parameters s x and s 3 control all high gradient phenomena such 

as shock waves and turbulence. These parameters as calculated from the changes of local 
Mach numbers and Reynolds (or Peclet) numbers within each element and are indicative 
of the actual local element flowfields. The contours of these parameters closely resemble 
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the flowfields themselves, with both and s 3 being large (close to unity) in regions of 
high gradients, but small (close to zero) in regions where the gradients are small. The 
basic role of s x and s 3 is to provide computational accuracy. 

(2) The second order implicitness parameters s 2 and s 4 are also flowfield dependent. 
However, their primary role is to provide adequate computational stability (artificial 
viscosity) as they were originally introduced into the second order time derivative terra of 
the Taylor series expansion of the conservation flow variables if . The primary role of s 2 
and s 4 is to provide computational stability. 

(3) The terms represent convection. This implies that if = 0 then the effect of convection 
is sraalL The computational scheme is automatically altered to take this effect into 
account, with the governing equations being predominantly parabolic-elliptic. Note that 
these effects are confined atl/* , not at U* . 

(4) The terms are associated with diffusion. Thus, with s 3 = 0, the effect of viscosity or 
diffusion is small and the computational scheme is automatically switched to that of Euler - 
equations where the governing equations are predominantly hyperbolic. 

(5) If the first order implicitness parameters s l and s 3 are nonzero, this indicates a typical 
situation for the mixed hyperbolic, parabolic and elliptic nature of the Navier-Stokes 
system of equations, with convection and diffusion being equally important. This is the 
case for incompressible flows at low speeds. The unique property of the FDMEI scheme 
is its capability to control pressure oscillations adequately without resorting to the 
separate hyperbolic elliptic pressure equation for pressure corrections. The capability of 
FDMEI scheme to handle incompressible flows is achieved by a delicate balance between 
s and .y, as determined by the local Mach numbers and Reynolds (or Peclet) numbers. If 
the flow is completely incompressible (M = 0), the criteria given by (19) leads to r, = l, 
whereas the implicitness parameter j 3 is to be determined according to the criteria given in 
(21). Make a note of the presence of the convection-diffusion interaction terms given by 
the product of bn in the s 2 terms and a b ; in the s 4 terms. These terms allow interactions 
between convection and diffusion in the viscous incompressible and/or viscous 
compressible flows. 

(6) If temperature gradients rather than velocity gradients dominate the flowfield, then s 3 is 
governed by the Peclet number rather than by the Reynolds number. Such cases arise in 
high speed, high temperature compressible flows close to the wall. 

(7) In the case of reacting flows the source terms B contains the reaction rates which are 
functions of the flowfield variables. With widely disparate time and length scales involved 
in the fast and slow chemical reaction rates of various chemical species as characterized by 
Damkohler numbers, the first order source term implicitness parameter s 5 is instrumental 
in dealing with the stiffness of the resulting equations to obtain convergence to accurate 
solutions. On the other hand, the second order source term implicitness parameter s 6 
contribute to the stability of solutions. It is seen that the criteria given by (27-28) will 
adjust the reaction rate terms in accordance with the ratio of the diffusion time to the 
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reaction time in finite rate chemistry so as to assure the accurate solutions with 
computational stability. 

(8) Various definitions of Peclet number and Damkohler numbers (Table 1) between the 
energy and species equations should be checked. V/hichever definition provides larger 
v alues of s 3 and s 5 must be used. The summary of the above definitions for implicitness 
parameters is shown in Table 2. 

(9) The tr ansitio n to turbulence is a natural flow process as the Reynolds number increases, 

causing the gr adien ts of any or all flow variables to increase. This phenomenon is the 
physical instability and is detected by the increase of s 3 if the flow is incompressible, but by 
both s 3 and fj if the flow is compressible. Such physical instability is likely to trigger the 
numerical instability, but will be countered by the second order implicitness parameters 
and/or s 4 to ensure numerical stability automatically. In this process, these flowfield 
dependent implicitness parameters are capable of capturing relaminarization, 
compressibility effect or dilatational turbulent energy dissipation, and turbulent unsteady 
fluctuations. > 

(10) An important contribution of the first order implicitness parameters is the fact that they can 
be used as error indicators for adaptive mesh generations. That is, the larger the 
implicitness parameters the higher the gradients of any flow variables. Whichever governs 
(largest first order implicitness parameters) will indicate the need for mesh refinements. In 
this case, all variables (density, velocity, pressure, temperature, species mass fraction) 
participate in resolving the adaptive mesh, contrary to the conventional definitions of the 
error indicators [10,15,16]. 

(11) Physically, the implicitness parameters will influence the magnitudes of Jacobians. Thus, 
Item(8) above may be modified so that the diffusion implicitness parameters s 3 and s 4 as 
calcuk^ from the Reynolds number and Peclet number can be applied to the Jacobians 
(a., b , c . corresponding to the momentum equations and energy equation, respectively. 
Furthermore, two different definitions of Peclet number (Pe,, Pe n ) would require the s 3 and 
s 4 as calculated from the energy and species equations to be applied to the corresponding 
terms of the Jacobians. Similar applications of the source term implicitness parameters s 5 
and s 6 should be followed for the source term Jacobian d with five different definitions of 
Damkohler number applied to the corresponding terms of d. In this way, high temperature 
gradients arising form the momentum and energy equations and the finite rate chemistry 
governed by the energy and species equations can be resolved accordingly. 

The FDMEI equations as given in (13) or (18) may be solved by either FDM or FEM. 
The standard linear Galerkin approximations of FEM lead to the results of central differences of 
FDM. However, the main difference between FDM and FEM arises when integration by parts is 
performed in FEM and the explicit terms of Neumann boundary conditions “naturally'’ appear as 
boundary integral forms. Thus, all Neumann boundary conditions can be directly specified at 
boundaries in FEM. This is not the case for FDM. Often, a rather cumbersome process must be 
taken for Neumann boundary conditions in FDM. 
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When dealing with all speed flow regimes such as in shock wave turbulent boundary layer 
interactions where compressible and incompressible flows, viscous and inviscid flows, and laminar 
and turbulent flows are intermingled throughout the flowfield domain, a computational scheme 
intended for only one type of flow physics and that does not account for other types of flow 
phenomena will fail. For example, the flow close to the wall in shock wave turbulent boundary 
layer interactions is incompressible (M < 0.1), whereas away from the wall the flow is 
compressible (supersonic or hypersonic). In this case, viscous flows change to inviscid flows. In 
between these two extremes the flowfield changes continuously, oscillating back and forth across 
the boundary layers of velocity and entropy, and leading edge and bow shocks. At any given 
computational nodal point or element, gradients of each variable (density, pressure, velocity, and 
temperature) may be very small or very large, so large that practically all currently available 
computational methods may fail. In order to succeed, it is necessary that the current flow physics 
everywhere be identified and so recognized, with specific computational schemes accorded to 
each and every computational nodal point and element. It is clear that such accommodations are 
available in (13) or (18). 

3. IMPLEMENTATION AND COMPUTATIONAL PROCESS 

As stated earlier, the governing equations for the Taylor series-modified Navier-Stokes system of 
equations, (13) or (18) may be applied to either FDM or FEM, or to the finite volume method 
(FVM). For application to FDM any of the existing finite difference methods can be used to 
obtain the standard finite difference analogs for either (13) or (18) such as the central schemes, 
upwinding schemes, or TVD schemes. The role of the FDMEI is them to enhance the 
computational accuracy above and beyond the limit of the current FDM capabilities. 

For applications to FEM we begin by expressing the conservation and flux variables and 
source terms as a linear combination of trial functions O a with the nodal values of these variables. 

u(i,<) = «„(i) u„« , r,H=*.WF.W 


G,(x,/)-4>.(i)G«(0 , B(*,i)-®.W B a (r) 

(29) 

Applying the generalized Galerkin approximations to (1 1) we obtain 


J<D a R(U,F„G„B)dQ = 0 

Q 

(30) 

or 



(31) 

where 



(32) 
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n„ = 8 „ - S i bld„ - s t i A t 2 d„d a 


(33) 
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or 




Oa <& p 


(36) 


+<*)*■ M"‘‘ ir 


where all Jacobians must be updated at each time step and 4> a represents the Neumann boundary 

trial and test functions, with a, (3 denoting the global node number and r, s providing the number 
of conservation variables at each node. For three dimensions, i, j = 1,2,3 associated with the 
Jacobians imply directional identification of each Jacobian matrix (a t , a^, a,, b p b 2 , b 3 , c u , c 12 , 
C| 3 , Cjj, c 22 , Cjj, Cjj, Cj 2 , C33) with r, s = 1, 2, 3, 4, 5 denoting entries of each of the 5 x 5 Jacobian 
matrices. 

It is important to realize that the integration by parts as applied to the generalized Galerkin 
approximations in FEM produces all Neumann boundary integrals. It is particularly advantageous 
that Neumann boundary conditions through re-evaluation of Jacobians normal to the boundary 
surfaces can simply be added to the boundary nodes for the stiffness matrix B in (34). On the 
other hand, all Neumann boundary conditions which appear in (36) act as source terms. These 
fe atures are absent in FDM, but implementations of Neumann boundary conditions can be handled 
by devising special forms of finite differences at boundary nodes. 

The gener alize d Galerkin approach of (13) may be replaced by the generalized Petrov- 
Galerkin methods. This process will require the RHS of (34) and (35) to be revised by adding, 
respectively. 
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with E being the quantities inside of the brackets of the convection terms in (34) and E„ are 
those in (34) and (3g are the Petrov-Galerkin parameters as defined in (A20), Appendix A. The 
role of the FDMEI is, as in FDM of different schemes, to enhance the solution accuracy above 
and beyond the Petrov-Galerkin methods. Otherwise, the formulation given by (13) represents 
the generalized Taylor-Galerkin method with accuracy enhanced by the FDMEI scheme. 


Similar results are obtained either by FDM or FEM with accuracy of computations derived 
primarily from the FDMEI equations of (13) or (18). However, with the increase of Reynolds 
number (say around Re »10 8 ), it is possible that accuracy may increase with applications of 
special functions such as Legendre polynomials of high degree modes characterizing extremely 
small turbulent microscales. Implementation of such high frequency modes can be achieved by ^ 
placing these modes between the comer nodes of isoparametric finite elements. Adaptively, such 
high modes can be chosen as needed for the resolution of turbulent microscales. Once again the 
diffusion implicitness parameter s 3 will play a crucial factor in determining the required degrees of 
Legendre polynomial The use of Legendre polynomial spectral modes superimposed onto 
isoparametric elements has been discussed in [9,12,15]. Its merit, however, has not yet been fully 
established for general applications. 


For turbulent flows with an extremely high Reynolds number, the phase error of the short 
wavelengths can be very large. In this case, it is necessary to add numerical dissipation terms to 
damp out the short wavelengths. Such numerical viscosities are conceptually different from the 
second order implicitness parameters whose role is to ensure stable solutions while preserving the 
solution accuracy dictated by the first order implicitness parameters. Toward this end, it is 
desirable to revise (18) in the form 




3E, 

dx 


dXjdXj y 


AU" +I 


= _Q-_Q- 


(37) 


where Q* is the numerical dissipation vector in terms of the second order tensor of numerical 
dissipation,^. , associated with the second order derivatives oftf. 


p a 2 ir _ A . a 2 u- 

= S :1 - — r — = p. Ax ; Ax ; ■ 


9 dx^Xj 


dxflxj 


(38) 


with [L being the numerical dissipation constant chosen as 0< p. < p o , where p a is set 
approximately equal to 2, but adjusted from numerical experiments. Note that the Galerkin 
integral of (38) (integration by parts once) leads to the first derivative of the trial and test 
functions combined with the nodal values of U£. In addition, note that the damping provided by 
the second order derivatives will not disrupt the formal accuracy of the FDMEI scheme. This 
process may be applied to (13) as well. 
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One of the most significant aspects of the FDMEI scheme is that for low Mach numbers 
(incompressible flow) the scheme will automatically adjust itself to prevent pressure oscillations. 
This adjustment is analogous to the pressure correction scheme employed for incompressible 
flows. Otherwise, the FDMEI scheme is capable of shock wave resolutions at high Mach 
numbers, and particularly well suited for dealing with interactions between shock waves and 
turbulent boundary layers where regions of high and low Mach numbers and Reynolds numbers 
coexist. In this case, the inviscid and viscous interactions are allowed to take place. To this end 
the second order implicitness parameters play the role of artificial viscosity needed for shock wave 
resolutions in the presence of flow diffusion due to physical viscosity. 

In order to understand how the FDMEI scheme handles computations involving both 
compressible and incompressible flows fundamental definitions of pressure must be recognized. 
Consider in the following that the fluid is a perfect gas and that the total energy is given by 

E = c r~£ + -U ( v,. (39) 

P 2 

The momentum equation for steady state incompressible rotational flow may be integrated to give 

J(p + 7P v , v ,),A = JK V Uf + i V ;.y<) + P E y* V ,®*]A 

P+iP v ; v ; =Po+ w (40) 

with 

w + -k,.) + p e -/» v /°*]a 

where co t is the component of a vorticity vector, p a is the constant of integration, and m denotes 
the spatial dimension. 

Combining (39) and (40) leads to the following relationship: 

p„=p(c,r + v,v,-£)-lV (41) 

tf Po as given by (41) remains a constant, equivalent to a stagnation pressure, then the 
compressible flow as assumed in the conservation form of the Navier-Stokes system of equations 
has now been turned into an incompressible flow, which is expected to occur when the flow 
velocity is sufficiently reduced (approximately 0.1 < M < 0.3 for air). Thus, (41) may serve as an 
equivalent equation of state for an incompressible flow. This can be identified element by element 
for the entire domain. Note that conservation of mass is achieved for incompressible flows with 
p o in (41) being constant, thus keeping the pressure from oscillating. 

Once the Navier-Stokes solution via FDMEI is carried out and all flow variables 
determined, then we compute fluctuations/ of any variable/. 


/' = /-/ 


( 42 ) 
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where / and / denote the Navier-Stokes solution and its time average, respectively. This 
process may be replaced by the fast Fourier transform of the Navier-Stokes solution. Unsteady 
turbulence statistics (turbulent kinetic energy, Reynolds stresses, and various energy spectra) can 
be calculated once the fluctuation quantities of all variables are determined. 

Before we demonstrate numerical examples, let us summarize why the FDMEI scheme is 
capable of handling low speed and high speed and compressible and incompressible flows, 
including shock waves and turbulent flows: (1) How is the transition from incompressible flow to 
compressible flow naturally and automatically accommodated without using two separate 
equations or two separate codes? This process is dictated by the first order convection 
implicitness parameter si as reflected by the Mach number changes and the expression of the 
stagnation pressure. (2) How is the shock wave captured? As the Mach number increases and 
its discontinuity is abrupt, the s 2 terms associated with second order derivatives together with 
squares of the convection Jacobian provide adequate numerical viscosities through second order 
derivatives, similarly as the Lax-Wendroff scheme. (3) How is the transition from laminar to 
turbulent flows naturally and automatically accommodated? This process is governed by the first 
and second order diffusion implicitness parameters (s 3 and sj as calculated from the changes of 
the Reynolds number. The terms associated with s 3 and s 4 are responsible for fluctuations of 
velocities, with the values of these implicitness parameters increasing with intensities of turbulence 
in conjunction with the diffusion gradient Jacobian and the squares of the diffusion Jacobian. This 
process allows the Navier-Stokes solutions to contain fluctuations which can be extracted by 
subtracting the time averages of the Navier-Stakes solutions. (4) How do the interactions 
between convection and diffusion take place? Changes of Mach numbers and Reynolds numbers 
as reflected by both convection and diffusion implicitness parameters close to the wall contribute 
to the uns teadin ess. Away from the walk they contribute to the transition between incompressible 
to compressible flows. (5) How are the stiff equations arising from widely disparate reaction 
rates of all chemical species treated? The most crucial aspect of the FDMEI scheme is its 
capability to identify the ratio of the resident time to the reaction time as calculated from five 
different definitions of the Damkohler numbers between the adjacent nodal points and time steps 
as reflected in the calculated first order implicitness parameter, s 5 , and the second order 
implicitness parameter, s 6 . These parameters provide precise degree of computational implicitness 
at every nodal point and every time step, contributing to the determination of accurate chemical 

reactions. 

4. APPLICATIONS 

We examine here various example problems:(a) flow over a flat plate, (b) shock wave 
turbulent boundary layer interactions on a compression comer, (c) 3D duct flows, and (d) lid- 
driven cavity flow. Linear isoparametric finite elements are used for the example problems. 

(a) Behavior of Flow field Dependent Implicitness Parameters on Flat Plate 

First of alL our concern is to test the behavior of FDMEI and FDMEI-FEM. Toward this 
objective we examine the flow over a flat plate investigated earlier by Carter [17] as shown in Fig. 
4a. The initial setting for the implicitness parameters are determined from the initial conditions of 
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the flowfield and subsequently updated after each time step until the steady state solution is 
reached. 

Corresponding to the mesh refinements and the flowfields at steady state shown in Fig. 4b, 
c, d, the contours of implicitness parameters s { and s 3 are given in Fig. 5. It is seen that the 
implicitness parameters themselves closely resemble the flowfield. There are little or no changes 
in Mach numbers and Reynolds numbers between adjacent nodes or elements far away from the 
surface of the plate as indicated by s, = r 3 = 0. ! Along the leading edge shock and boundary 
layers, both and s 3 move toward unity indicating that gradients of all variables increase. The 
final flowfields, as shown in Fig. 4b,c,d, are the consequence of these implicitness parameters. 
The implicitness parameters s 2 and s 4 are the compliances of and s v respectively, with their 
primary roles being the artificial viscosity. Thus, the first order implicitness parameters (s v s 2 ) 
help to resolve the high gradients ensuring the accuracy of the solution. While on the other hand, 
the second order implicitness parameters (s 2 , s A ) ensure computational stability. 

Compu tatio ns of wall pressure, wall skin friction, u-velocity, v-velocity, density and 
temperature distribution are shown in Figs. 6a through 6f. The comparison with the Carter's data 
indicates reasonable agreements. 

(b) Supersonic Flow on a Compression Comer 

In this example we demonstrate calculations of supersonic flow on a compression comer. 
The inlet boundary conditions (non-dimensionalized) are p = 1, M = 2.25, p = 0.14, Re = 105, 
Pr = 0.72, and v = 0, with adiabatic wall condition. The steady state background mean flowfields 
for the compression comer are shown in Fig 7a. In these calculations, all perturbation 
(fluctuation) variables are determined from time averages of the Navier-Stokes solutions 
according to (35). The horizontal and vertical perturbation velocities (u, v') at locations close to 
the wall (jc = 0.10256 m, y = 0.001 m) and away from the wall (x = 0.10256 m, y - 0.04 m) are 
shown in Fig. 7b. Note that u is extremely unsteady whereas v' is significantly less unsteady 
close to the walL Away from the wall, both u and v' are almost steady. These trends are 
reflected in the turbulence (Reynolds) stresses as shown in Fig 7c. Turbulent kinetic energy 
distributions at the locations upstream of the comer (x = 0.0513 m) and downstream of the comer 
(x = 0.1333 ra) are shown in Fig 7d. We observe that the turbulent kinetic energy downstream of 
the comer is s i gnific antly larger than the upstream. No turbulent statistics calculations (wave 
numbers or frequencies vs power spectral density) are attempted at this time as turbulence 
microscales are not resolved in this example. 

It should be noted that the above results obtained without turbulence models or without 
the standard DNS solutions (neither spectral nor DNS-mesh refinements) are regarded as the 
consequence of the time- averaging of the FDMEI Navier-Stokes solutions. This implies that the 
fluctuation of variables between nodal points (Fig 2a) and between time steps (Fig 2b) as reflected 
in terms of the implicitness parameters (j.) have contributed to these physical phenomena, with 
compressibility and shock waves dictated by the Mach number-dependent s v and with 
incompressibility and turbulent fluctuations dictate by the Reynolds number or Peclet number- 
dependent (s 3 ). An equal participation of r, and r 3 will be responsible for shock wave turbulent 
boundary layer interactions. A comparison of the results of the FDMEI scheme with the K-e 
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turbulent model and experimental data is shown in Fig 7e. It is seen that the FDMEI results 
compare more favorably with those of measurements [18]. 

(c) FDMEI Analysis of Three Dimensional Flows 

To demonstrate the effectiveness of the flowfield-dependent implicitness parameters in 3- 
D flows at the steady state, we examine the spatially evolving boundary layer (Figs. 8a through 
8e). Note that the contours of s x and s 3 (Fig. 8c) show the boundary layer effects in which both s x 
and Sj are indicative of rapid changes of Mach numbers and Reynolds numbers respectively, larger 
(close to unity) on the wall, but small (closer to zero) away from the walL The velocity vectors 
and RMS error distributions versus interactions are shown in Figs. 8d and 8e, respectively. 

(d) Demonstration of Compressibility vs Incompressibility 

We ask the question: Can a single formulation or computer program originally designed 
for high speed compressible flows be applied to analyze the low speed incompressible flows? The 
advantage of FDMEI is to respond positively to this question. To prove the point, let us examine 
the lid-driven cavity flow at the steady state (Figs. 9a through 9f). Notice that, for M = 0.1, 
density changes occur closer to the lid, whereas, for M = 0.01, density is constant throughout the 
domain (Fig. 9e), corresponding to P„ being variable and constant, respectively (see Eq. (32)). 
The equation of state for compressible flows is automatically switched over to accommodate the 
incompressible flows. This advantage is contrary to the previous practice such as the Table look- 
up for the equation of state for incompressible flow handled separately through hyperbolic elliptic 
equation as derived from the continuity equation combined with the momentum and energy 
equations. Comparisons of the results of FDMEI with those of the independent incompressible 
flow code of Ghia et al [19] are very favorable as shown in Figs. 9a through 9f. 

5. CONCLUDING REMARKS 

The validity of the proposed new approach to computational fluid dynamics has been 
demonstrated through some example problems. Excluded from these examples are reacting flows 
which are reported elsewhere [16]. Also excluded is the effect of additional spectral modes of 
Legendre polynomials which are described in [15]. None of the example problems have been 
carried out with mesh refinements required for resolving turbulent microscales due to the 
limitation of computer time. The following concluding remarks are provided: 

(a) The flowfield-dependent implicitness parameters as calculated from the current 
flowfield information are indicative of the magnitude of gradients of all variables and 
adjust the computational schemes accordingly for every nodal point or element, rather 
than dictated by arbitrarily selected constant parameters throughout the domain. 

(b) The first order implicitness parameters s v s v and s s as calculated from Mach numbers, 
Reynolds or Peclet numbers, and Damkohler numbers, respectively, ensure the 
solution accuracy, whereas the second order implicitness parameters s 2 , s 4 , and s 6 
which are determined as compliances of s v s v and s y respectively, assist in the 
solution stability. 

(c) The FDMEI method is capable of resolving mutual interactions and transition between 



19 


viscous and inviscid flows, compressible and incompressible flows, and laminar and 
turbulent flows, in all speed regimes. 

(d) Further research on FDMEI is required in order to investigate many other physical 
phenomena including hypersonics and reacting flows with high temperatures in 3D 
geometries. 
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Table 1 Definitions of Nondimensional Flowfield Quantities 


p(v . V)v = -Vp + p[V 2 v +^V(V • v)] 
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r * , 
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Mach number 

M 
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a 

A inertial force 
B pressure force 

Reynolds number 

Re 

p uL 

A inertial force 
C viscous force 

Peclet number, I 

Pe, 

putept 

k 

E convective heat transfer 
F conductive heat transfer 

Peclet number, II 

Pen 

uL 

D 

I convective mass transfer 
J diffusive mass transfer 

Damkohler number, I 

Da, 

Lw t 

K mass source 
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I convective mass transfer 

Damkohler number, II 

Da„ 
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K mass source 

J diffusive mass transfer 

Damkohler number, HI 

Da ra 

qL 

N heat source 

Hu 

E convective heat transfer 

Damkohler number, IV 

Da iv 

2! 

kT 

N heat source 

F conductive heat transfer 

Damkohler number, V 

Da v 

qL l 
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N heat source 

G diffusive heat transfer 



























Table 2 


Flowfield Dependent Implicitness Parameters 


First order convection implicitness parameter 


min(r,l) 

0 
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r> a 

r < a , *0 

M min =0 


r = J M 


-mL/m, 


Ensures solution accuracy. 
Strongly flowfield dependent, with 
high gradients characterized by large 
Mach number changes between 
nodal points or within element and 
between time steps. 


- Second order convection implicitness parameter 



Ensures solution stability. 
Flowfield dependent artificial 
viscosity for convection process 


s 3 - First order diffusion implicitness parameter 



Ensures solution accuracy. 

Strongly flowfield dependent, with 
high gradients characterized by large 
changes in Reynolds number or 
Peclet number between nodal points 
or within element and between time 
steps. Diffusion gradient behavior 
may be dictated by Peclet number 
when temperature gradients are 
high. Choose whichever (Reynolds 
or Peclet number) provides the largo* 
value for s* 


s 4 - Second order diffusion implicitness parameter 


: . 0 < n < 1 


Ens ures solution stability. 
Flowfield dependent artificial 
viscositv for diffusion process 


s 5 - First order source term implicitness parameter 

Ensures solution accuracy. 
Strongly flowfield dependent, with 
high reaction rate gradients 
characterized by large Damkohler 
number changes between nodal 
points or within element and 
between time steps. 
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s - Second order source term implicitness parameter 
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" , 0 < n < 1 


Ens ures solution stability. 
Flowfield dependent artificial 
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Fig 1 Supersonic and Hypersonic Flows 
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Fig 2 Spatial and T emp oral Flowfidd Dependent Implicitness Parameters 
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fig 3. Relationships between the first and second order 

implicitness parameters. Stable solutions occur in the 
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second order implicitness parameters to preserve the 
solution accuracy as dictated by the first order 
implicitness parameters. 
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Fig. 4 Flat plate problem - initial and adaptive meshes and their corresponding density 
con t o urs 









1 — ^ 

(a) s, contours 



(b) Sj contours 


Fig. 5 Flown eid-depenrient first order convection and diffusion implicitness 
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(b) Fluctuation velocities 
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fig. 7 Continued 
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APPENDIX A 

ANALOGIES OF FDMEI TO CURRENTLY AVAILABLE FDM AND FEM SCHEMES 

Analogies of FDMEI to currently available computational schemes of FDM and FEM are 
summarized below. 

A1. Analogies of FDMEI to FDM 


Some of the FDM schemes are compared with FDMEI in the following Table. 
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Other schemes of FDM are compared with FDMEI as follows: 

(a) Lax- Wendrojf Scheme 

The Lax-Wendroff scheme without artificial viscosity takes the form 



This scheme arises if we set in FDMEI 

a..=a._ l =a , j,= 0 , s 2 =0 , j 3 =0 , j 4 =0 

“ 2 


(b) Lax-Wendroff Scheme with Viscosity 

The Lax-Wendroff scheme with viscosity is given by 
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F ; , i + F A t 


_ Vi+L 


-^ a H (F«- F ,)+^(U M -0,) 


F + F,_, _Ar_ 
5 2 2Ajc 


f; , = i^L-^a, H (F ; -F m ) + DA U, - U M ) 


This scheme arises if we set 

D. , = D._l = as, , s 2 = 0 , S3 =0 , s 4 = 0 

,+ T J 2 

This implies that as, in FDMEI plays a role of artificial viscosity which is manually implemented in 
the Lax-Wendroff scheme. 

(c) Explicit McCormack Scheme 

Combining the predictor corrector steps of McCormack scheme we write 



The FDMEI becomes identical to this scheme with the following adjustments: 
a * = a 4 = a 

F," - F", = f ; +1 - F," + F.^ - F h 

s, = 0 , s 2 = 0 , S3 = 0 , s 4 =0 

and the s 2 term in the FDMEI method is equivalent to 

D. = |(U' H - 4U-„ +6U" - 4U"_, + 0;_,) 

This again is a manifestation that shows the equivalent of the s 2 terms is manually supplied in the 
McCormack method. 


(d) First Order Upwind Scheme 
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This scheme is written as 


Aur 1 = 


-fK- F 4) 

-^fi( F '" +F "' ) 4W u ''- u ')] 

- |(F'+F-,)-i|ai(u:-U' 1 )]} 


(A4) 


The FDMEI analogy is obtained by setting 

F’ = ^F- , , F-,=iF-, 

s 2 ac(Aur‘ - 2 au;;; + au-;‘) = | a|(u; +1 - u;. t ) 

where C is the Courant number. 

(e) Implicit McCormack Scheme 

With all second order derivatives removed from (11) we obtain the implicit McCormack 
Scheme by setting s t = l, s 2 = 0, s 3 = 0, s 4 = 0. However, it is necessary to divide the process into 
the predictor and corrector steps. Once again the flowfield-dependent implicitness parameters for 
FDMEI will allow the computation to be performed in a single step. 

(f) PISO and SIMPLE 

The basic idea of PISO and SIMPLE is analogous to FDMEI-FEM in that the pressure 
correction process is a separate step in PISO or SIMPLE, whereas the concept of pressure 
correction is implicitly embedded in FDMEI-FEM by updating the implicitness parameters based 
on the upstream and downstream Mach numbers and Reynolds numbers within an element. 

The elliptic nature of the pressure Poisson equation in the pressure correction process 
resembles the terms embedded in the B^ r , terms in (28a). Specifically, examine the s 2 terms 
involving a irq a m and b irq a M and s 4 term involving a irq b Pr All of these terms are multiplied by 
<h j which provide dissipation against any pressure oscillations. Question: Exactly when is such 
dissipation action needed? This is where the importance of implicitness parameters based on 
flowfield parameters comes in. As the Mach number becomes very small (incompressibility 
effects dominate) the implicitness parameters s 2 and s 4 calculated from the current flowfield will 
be i ndicativ e of pressure correction required. Notice that a delicate balance between Mach 
number (s 2 is Mach number dependent) and Reynolds number or Peclet number (s 4 is Reynolds 
number or Peclet number dependent) is a crucial factor in achieving a convergent and stable 
solutions. Of course, on the other hand, high Mach number flows are also dependent on these 
implicitness parameters. In this case all implicitness parameters, s t , s 2 , r 3 , *4 will play important 

roles. 
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A2. Analogies of FDMEI to FEM 


(a) Generalized Taylor Galerkin ( GTG ) with Convection and Diffusion Jacobians 

Earlier developments for the solution of Navier-Stokes system of equations were based on 
GTG without using the implicitness parameters. They can be shown to be special cases of 
FDMEI- FEM. 

In terms of both diffusion Jacobian and diffusion gradient Jacobian, we write 


dG, _ U 3U 3V, 
dt b ‘ dt Cii dt 


with 


_dG ± = dG, = dV 

bi " au ’ 9 av ; ’ ' dxj 


Thus it follows from ( 10) or ( 1 1) with Sj — — s 4 — — s 6 — 0 and — 1 that 


AU” 1 = a/-| 3— r^+B 
l dx, ox, 


V At 2 d f 3F ; dG l+B y 


2 dt 


V 


dx, dx, 


0(A,>) 


(A5) 


Using the de finitio ns of convection, diffusion, and diffusion rate Jacobians discussed in Section 2, 
the temporal rates of change of convection and diffusion variables may be written as follows: 


3F- _j 

\ 3LY - 


'JUL- 

dGj 

dt \ 

l ' at J 

a t 


dxj 


B 


n 





, d¥ i dG? „„ 
- U ) — - 1 T J —+ B" 

dx j ux i 


/j 


(A6) 



aGf 1 _ 

fb -*i] 

au *** _a_ 


dt 

I* 

At + a Xj 

l 4 At J 


Substituting (A6) and (A7) into (A5) yields 


(A7) 
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AU^ 1 =a/-|^ l --^ l +B 

V dx i J 

\ V WT rt+1 


-5, 


AU 

A/ 


Assuming that 


At\ 

a " 

2 

dx i 

B*+i 


dt 


e : = 

b- 




' sau- 1 3f; scr . 

—a t o 

V ' dx j dx J dx J ) 


*L. fl 


(A8) 


and neglecting the spatial and temporal derivatives of B, we rewrite (A8) in the form 


1- 


A/ 2 a 


2 dx, 


Al-Ll 


, a a • h — - I— — >AU = H 

l ; Aija^. 


(A9) 


( aF r . 

ac ( , b' 

f At 2 a 

' aF."| 

a ' 

1 a* 

ar, +B 

) 2 a*,. 

I'^J 


Here the second derivatives of G ; are neglected and all Jacobians are assumed to remain constant 
within an incremental time step but updated at subsequent time steps. 

Applying the Galerkin finite element formulation, we have an implicit scheme, 

(A*5„ + B^AlC = Hi + fC + Ni (A10) 


where 


B <*0« ~ -I J 


a L 


^ if ft 

a. a H — ■ — 

,r<r 131 At 


**a*9.i 


dQ. 


HI = ArJ +Gp j ,)-d> a ^>pBp r -^-a i „<I> 0 , j <I> pj F^j^ 

Q L J 


N" +l = A£lff a . a- +^- 
or 2 J[ ,rq ,sq At 


<b« Au:t‘n ; dr 


S.J t 


Ni =-J [Al<i>a(F,; + G' ) - 6. F'.,l n,dT 

r L -* 


Here we note that the algorithm given by (A10) results from (29) by setting 
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- _ c _ _n c - 1 h a = C. JAt, and neglecting the terms with b and derivatives of G, and 
B, the form identical to that reported by Hassan et al [13]. 

(b) GTG With Convection Jacobians 

Diffusion Jacobians may be neglected if their influences is negligible. In this case the 
Taylor-Galerkin finite element analog may be derived using only the convective Jacobian from the 
Taylor series expansion. 


v „ dU Ar 2 9 U 
U = U +Ar— — +— — r- 
dt 2 dt 


r + °M 


(All) 


where 


3U = J£._ B = -a, B 

dt dx, dx, ' dx, dx. 


d 2 U _ _d_f 3U + dG 
dt 2 dt l ' dx, dx, 


H 


(A12) 


or 


a 2 u 


dv) 

d 


hdxj 

1 

dx, 

l ' * »J 


d , DX dB 


(A13) 


Substituting (A12) and (A13) into (Al 1), we obtain 


AlT*=Ar 


dF f 


d 

( dV) 

d 2 (a i G i ) 

d 

L 

dr, 

■is: +B+ 2 


r^,J 

a xfr 

dx, 


dB 

dr 


(A14a) 


Expanding dFVdr at (n+1) time step 


dF,"* 1 _ 
dr 


dF, dG 


+ B 


J dAU"* 1 dF" dGf 1 , „„ 
' [ J dx J dx } dXj 


dxj dXj 

and substituting the above into (Al 1 - A13), we arrive at AU" +l in a form different from (A14a), 
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(d) Characteristic-Based Zienkiewics/Codina Scheme 

This scheme arises from Eq. (10) by splitting the FDMEI Navier-Stokes system of 
equations into three parts for continuity, momentum, and energy, separately. 


Continuity 


AU" +l = -At 


9f; 

3a ; AU" +1 

A t_ 

f 3 2 a f F" a 2 a,a -AU" +1 ^ 

3 — + Si 
a*, 

dx. 

2 

dxjdxj 2 Bx i dx j ^ 


(All) 


where all diffusion terms are neglected. Setting 


AU" +l Ap" +l 

f" -+ pv; u* 

r 1 a i AU" +l s,Apv” — > 0 ,Af 7 " 
■2 a fF" e lP "8 l7 
x, 2aj a ; .AU- +1 O^A/T'S* 


These substitutions to (A21) lead to 




Ap"" 1 =lp-ApJ = "At 


3 (p v ,)" ^(Apy,) 1 




3*; 




-Are, 


^2 n 

° p 

dXidx; 


+ e 2 


3 2 A p 
dx,dx 


#1+1 ^ 


. ) 


(A22) 


which is identical to (33) with (Apv,)^ 1 = AU; being the intermediate step in [14], This 
represents the pressure correction equation. 


Momentum 

( A P V ,) 




= -At 


a(p v ,v.)" ax; a p - 9Ap- +1 

+0 2 


3x, 


dx, dxj 


dx ; 


-(i-e,)[T-v,^(pv,v j + psj] 


(A23) 


which is similar to (30) of [14] with a, = v ;> 1 - 0 2 = s,, and all terms of r 2 , b,, and c.. being 
neglected in FDMEI. 
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Energy 


(ApE)" + ‘ = -A t 


3(pEv; + VjpY_ d 


dX; 


dX: 


, dT 

k- — +x„v ,• 

OX; 


IJ J 


(A24) 


_ Ar_ v <L 

2 ' dxjdx j 


(pEV; + VjP)" 


which is identical to (40) of [14] with alln + 1 terms being neglected in FDMEI. 

The solution steps begin with (A-23), foUowed by (A-21) and (A-24). Note that the 
pressure corrections for incompressible flows are internally carried out in FDMEI as the pressure 
second derivatives arise automatically in Ecj(lO). Note also that in FDMEI all implicit terms may 
be recovered if so desired. 
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AU 3AU A t 3 2 AU _ 

ha. a a = 0 

A t ' dx t 2 ‘ J dx$xj 


For the steady state non-increraental form in 1-D we write (A16) in the form 

du . o 2 d 2 u . 

a At r = 0 

dx 2 dx 2 

Taking the Galerkin integral of (A 17) leads to 


2 dx 2 


(A16) 


(A 17) 


or 


\wjf ) a^-dx = 0 

J OX 

for vanishing Neumann boundaries. Here wjf' 1 is the Petrov-Galerkin test function, 

Wj‘ ) =oy+ ah^- 


(A18) 


(A19) 


with a = C/2 and C = aAt/Ax being the Courant number. 

For isoparametric coordinates in two dimensions, the Petrov-Galerkin test function 
assumes the form 




af'i 


(«) 


(A20) 


with being the Petrov-Galerkin parameters 


a? = cothf-^- 


- a n = coth] 
"5 


ffl-i 


8i = 


V, 


where R is the Reynolds number or Peclet number in the direction of isoparametric coordinates 
(£, ti). Note that the GPG process given by (A16)-(A20) leads to the Streamline Upwinding 
Petrov-Galerkin (SUPG) scheme as a special case. 
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AU** 1 = Ar 


( 


0F, dG : 


v 3*i 


+ B 


+ 


Ar 2 

\±( 

2 * 

[at,. 


dAU^ 3F- 
a a t-a- — — 


M. 4 M r,» 


obc.-obc^ cte; 


dr 


H" = 


1- 


Ar 2 a 


2 dX; 


a 


-iMj 


AU" 


H" = Ar 


f ». 


fB T + Ax‘ a 

f aF.^1 

A ' 

< dx, 

dX; 

J 2 ax, 



(A 14b) 


(A14c) 


where second derivatives of G, is assumed to be negligible and B is constant in space and time, 
arriving at an implicit finite element scheme, 

(a*s„ + b^„)au;;‘ = h-„ + n"' + n;, (Ai5) 


where 


A*=j<wn 

a 




d> d> a 

'•'a i^P,/ 


d£l 


h: 






+ G 


a L 




d) d>» Fa" 
irs^ a ,1^ 3 , P jt 


dCl 



n,dT 


It should be noted that the form (A14c) arises from (25) with s, = s 3 = s 4 = b- = 0 and s 2 = 1, an 
algorithm similar toHassan et al [13]. 


(c) Generalized Petrov-Galerkin (GPG) 

The Generalized Petrov-Galerkin (GPG) method can be identified by setting $j = s 2 = 1, 
s i = s 4 = 0, b ; = c.. = d = 0, Q" = 0, E ; = a,, and E, y = jAr 2 a ; a y , so that (11) takes the form 



